Non-stationary heat conduction in one-dimensional chains with conserved momentum. 
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The Letter addresses the relationship between hyperbolic equations of heat conduction and mi- 
croscopic models of dielectrics. Effects of the non-stationary heat conduction are investigated in two 
one-dimensional models with conserved momentum; Fermi-Pasta-Ulam (FPU) chain and chain of 
rotators (CR). These models belong to different universality classes with respect to stationary heat 
conduction. Direct numeric simulations reveal in both models a crossover from oscillatory decay of 
short-wave perturbations of the temperature field to smooth diffusive decay of the long-wave per- 
turbations. Such behavior is inconsistent with parabolic Fourier equation of the heat conduction. 
The crossover wavelength decreases with increase of average temperature in both models. For the 
FPU model the lowest order hyperbolic Cattaneo-Vernotte equation for the non-stationary heat 
conduction is not applicable, since no unique relaxation time can be determined. 



It is well-known that parabolic Fourier equation of heat 
conduction implies infinite speed of the signal propaga- 
tion and thus is inconsistent with causality [1, 2, 3, 4, 5]. 
Numerous modifications were suggested to recover the 
hyperbolic character of the heat transport equation [2]. 
Perhaps, the most known is the lowest-order approxima- 
tion known as Cattaneo-Vernotte (CV) law [1, 2]. In its 
one-dimensional version it is written as 

(1 + T^)g=-A^VT (1) 

where k is standard heat conduction coefRcient and r 
is characteristic relaxation time of the system. The lat- 
ter can be of macroscopic order [5]. Importance of the 
hyperbolic heat conduction models for description of a 
nanoscale heat transfer has been recognized [6, 7]. 

Only few papers dealt with numeric verification of such 
laws from the first principles [8] . As it is well-known now 
from numerous numeric simulations and few analytic re- 
sults, the relationship between the microscopic structure 
and applicability of the Fourier law for description of the 
stationary heat conduction is highly nontrivial and de- 
pends both on size and dimensionality of the model [9]. 
In particular, the heat conduction coefficient can diverge 
in the thermodynamic limit. Hyperbolic equations de- 
scribing the non-stationary heat conduction inevitably 
include more empiric constants and therefore their re- 
lationship to the microscopic models can be even less 
trivial. To the best of our knowledge, no conclusive data 
exist in this respect. 

This Letter deals with a study of spatial and tempo- 
ral peculiarities of the non-stationary heat conduction in 
two simple one-dimensional models with conserved mo- 
mentum - Fermi-Pasta-Ulam (FPU) chain and chain of 
rotators (CR). From the viewpoint of the stationary heat 
conduction, these two systems are known to belong to 
different universality classes. Namely, in the FPU chain 
the heat conduction coefficient diverges with the size of 
the system [10], whereas in the CR model it converges to 
a finite value [11, 12, 13]. So, it is interesting to check 



whether other differences between these models models 
will reveal themselves in the problem of non-stationary 
heat conduction. 

In order to investigate this process, one should choose 
the parameters to measure. This question is not easy, 
since the situation in this problem is different from the 
stationary heat conduction, where only one commonly 
accepted macroscopic equation exists and only one em- 
piric parameter should be computed. The simplest CV 
law already has two independent coefficients, whereas 
more elaborate approximations can include even more 
parameters. Just because many different empiric equa- 
tions exist, it is not desirable to pick one of them ab ini- 
tio and to fit the data to find particular set of constants. 
Instead, it seems reasonable to look for some quantity 
which will characterize the process of the non-stationary 
conduction and can be measured from the simulations 
without relying on particular approximate equation. For 
this sake, we choose the characteristic length which char- 
acterizes the scale at which the nonstationarity effects 
are significant. 

In order to explain the appearance of this scale, let us 
refer to ID version of the CV equation for the tempera- 
ture: 

where a is the temperature conduction cocfRcicnt. 

Let us consider the problem of non-stationary heat 
conduction in a one-dimensional specimen with periodic 
boundary conditions T{L,t) = T{0,t), where T{x,t) is 
the temperature distribution, L is the length of the spec- 
imen, t > 0. If it is the case, one can expand the tem- 
perature distribution to Fourier series: 

oo 

T{x,t)^ an{t) eyi-p{2iTinx / L) (3) 

n— — oo 

with an{t) = al„(t), since T{x,t) is real function. 
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Substituting (2) to (3), one obtains the equations for 
time evolution of the modal amplitudes: 

rdn + a,i + in^n^aan/ L'^ = 0. (4) 
Solutions of Eq. (4) are written as: 

a„(i) = Cin exp(Aif) + C2„(A2t) 

Ai,2 = (-1 ±x/l~ WiT^n^aT/L^^ /2 

where Ci„ and C2n are constants determined by the ini- 
tial distribution. 

/^From (5) it immediately follows that for sufhciently 
short modes the temperature profile will relax in oscilla- 
tory manner: 

CLn{t) exp(— t/2T) exp(ia;„t), (6) 

If the specimen is rather long (L >> ATTy/ar) then for 
small wavenumbers (acoustic modes): 

Ai w -1/r, A2 ~ -An^n^a/L^. (7) 

The first eigenvalue describes fast initial transient re- 
laxation, and the second one corresponds to stationary 
slow diffusion and, quite naturally, does not depend on 
T. So, we can conclude that there exists a critical length 
of the mode 

r = 47rV^, (8) 

which separates between two different types of the relax- 
ation: oscillatory and diffusive. The oscillatory behavior 
is naturally related to the hyperbolicity of the system. 
Existence of this critical scale characterizes the deviance 
of the system from parabolic Fourier law. 

This critical wavelength scale I* can be measured di- 
rectly from the numeric simulation without relying on 
any particular empiric equation of the non-stationary 
heat conduction. The numeric experiment should be de- 
signed in order to simulate the relaxation of thermal pro- 
file to its equilibrium value for different spatial modes of 
the initial temperature distribution. We simulate the pe- 
riodic chain of particles with conserved momentum with 
Hamiltonian 

^ 1 

n=l 

In order to obtain the initial nonequilibrium temperature 
distribution, all particles in the chain were embedded in 
the Langevin thermostat. For this sake, the following 
system of equations was simulated: 

iln = V' {Un+l - Un) - V {Un - Un-l) - ^nU„ + in 

n = l,...,iV (10) 




Figure 1: Relaxation of initial periodic thermal profile in the 
chain of rotators, Z ^ 64, L = 1024, (a) To = 0.2, A = 0.05 
(oscillatory decay) and (b) To — 0.5, A = 0.15 (smooth decay 
of the initial thermal profile). 

where 7„ is the relaxation coefficient of the n-th particle 
and the white noise ^„ is normalized by the following 
conditions: 

iCn) = 0, (61(^1)6(^2)) = 2/^r„<5„fc5(tl - t2), (11) 

where T„ is the prescribed temperature of the n-th par- 
ticle. The numeric integration has been performed for 
7„ = 0.1 for every n and within time interval t = 250. 
After that, the Langevin thermostat was switched off and 
relaxation of the system to a stationary temperature pro- 
file was studied for various initial distributions T„ for two 
particular choices of the nearest-neighbor interaction de- 
scribed above (FPU and chain of rotators): 

Vi{x) = x^/2 + x^/4:, ^2(2;) = l-cosx. 

Separate analysis of individual spatial modes will pro- 
vide insight into the global behavior of the system only if 
these modes are, at least approximately, not interacting. 
For CV equation (2) this is the case since it is linear. 
However, we do not rely on it a priori and the absence 
of interaction between different spatial relaxation modes 
should be checked numerically. We simulate the relax- 
ation of the initial thermal profile comprising five differ- 
ent modes: 

8 

Tn=To + Y, A, cos[27r(n - l)/2'] (12) 

i=4 



Figure 2: (Color online) Evolution of the relaxation profile in 
the chain of rotators with change of the mode length Z. Time 
dependence of the mode maximum r(l + Z/2) (red lines) and 
minimum r(l) (blue lines) are depicted with average temper- 
ature To = 0.4 and Z = 16 x 2''~^, k = 1, 5, scaling time 
tk = 2'°~^. For all simulations length of chain L — 1024. 

in cyclic chain of = 2^^ particles, with average temper- 
ature To = 1 and modal amplitudes A4 ~ ... ~ = 0.02 
for the CR potential, and with Tq ~ 20, A4 = ... = Ag = 
0.4 for the FPU potential. Both simulations demon- 
strated almost complete lack of interaction between the 
modes. No other modes were excited with visible am- 
plitudes. Each mode in the collective excitation relaxed 
similarly to the profile obtained when it was excited indi- 
vidually. So, it is possible to conclude that for given val- 
ues of the parameters the equation of the non-stationary 
heat conduction should be approximately linear and sep- 
arate analysis of spatial relaxation modes is justified. 

In order to study the relaxation of different spatial 
modes of the initial temperature distribution, its profile 
has been prescribed as 

T„ ^Ta + A cos[27r(7i -1)/Z] (13) 

where Tq is he average temperature, A - amplitude of 
the perturbation, Z - the length of the mode (number of 
particles). The overall length of the chain L has to be 
multiple of Z in order to ensure the periodic boundary 
conditions. The results were averaged over 10^ realiza- 
tions of the initial profile in order to reduce the effect of 
fluctuations. 

Typical result of the simulation is presented at Fig. 1. 
The chain of rotators of the same length N = 1024 and 
the same modal wavelength Z = 64 demonstrates qual- 
itatively different relaxation behavior for different tem- 
peratures - the oscillatory one for lower temperature and 
the smooth decay - for higher temperature. This obser- 
vation suggests that the critical wavelength mentioned 
above, if it exists, should decrease with the temperature 



Figure 3: (Color online) Evolution of the relaxation profile 
in the FPU chain with change of the mode length Z. Time 
dependence of the mode maximum T(l + Z/2) (red lines) and 
minimum r(l) (blue lines) are depicted with average temper- 
ature To = 20 and Z = 16 x 2''^-^, k = 1, 7, scaling time 
tk = 2''-\ length L = 1024. 

increase. However, its existence should be checked for 
constant temperature and varying wavelength. 

Such simulations are presented at Fig. 2 (for the CR) 
and Fig. 3 (for the FPU chain). In both models one 
observes oscillatory decay for the short wavelengths and 
smooth exponential decay for relatively long waves. It 
means that for both models there exists some critical 
wavelength I* which separates two types of the decay and 
thus the effect of the non-stationary heat conduction is 
revealed. 

Results presented at Fig. 3 allow one to conclude that 
the critical wavelength for the FPU chain for given tem- 
perature may be estimated as 512 < /* < 1024. The 
interpretation of Fig. 2 is not that straightforward. It is 
clear that 32 < T < 128, but for Z = 64 the result is not 
clear. Within the accuracy of the simulation, it seems 
that only finite number of the oscillations is observed. It 
is possible to speculate that such behavior is not consis- 
tent with the lowest order CV equation, since expressions 
(6), (7) suggest either infinite number of the oscillations, 
or at most single crossing of the average temperature or 
no crossing at all. Possible interpretation may be that if 
the modal wavelength is close enough to the critical, the 
second-order CV model is not sufficient any more and 
the nonlocal effects of higher order should be taken into 
account. Still, these conclusions should be verified by 
more detailed simulations in the vicinity of the crossover 
wavelength. 

The latter observation has motivated us to check 
whether the data of numeric simulations in these one- 
dimensional models offer a support for the CV macro- 
scopic equation. For this sake, one can check another 
prediction of this equation - the independence of the am- 
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Figure 4: (Color online) Exponential decay of the normal- 
ized oscillation amplitude in the chain of rotators A{t) — 
{T{1 + Z/2)(t) -To) /A) for the average temperature To = 0.3, 
initial amplitude A = 0.05 and different periods of the ther- 
mal profile Z = 16 X 2*"^ k = 1,2,3. The straight lines 
illustrates the decay of the maximum envelope according to 
A = exp( — At) with universal value A — 0.015 for all three 
simulations. 

0. 




Figure 5: (Color online) Exponential decay of the normalized 
oscillation amplitude in the FPU chain a{t) = {T{l+Z/2){t)- 
To)/A for the average temperature To = 10, initial amplitude 
A = 0.05 and different periods of the thermal profile Z = 
16 X 2*^^^, k — 1, 5. The straight lines illustrates the decay 
of the maximum envelope according to A = exp(— At) with 
values A = 0.0015, 0.003, 0.008, 0.024 and 0.06 for fc = 1, 2, 
3, 4 and 5. 

plitude decrement of the relaxation profile on the wave- 
length in the oscillatory regime (6). The results of sim- 
ulation are presented at Fig. 4 (CR) and Fig. 5 (FPU). 
One can see that for the chain of rotators the above pre- 
diction more or less corresponds to the simulation results. 
For the FPU chain the decrement is strongly dependent 



on the wavelength, at odds with the CV equation. In 
this latter case, no unique relaxation time exists. 

To summarize, we reveal the hypcrbolicity effects of the 
non-stationary heat conduction in one-dimensional mod- 
els of dielectrics without relying on any particular em- 
piric equation. There exists a critical modal wavelength 
/* which separates between oscillating and diffusive re- 
laxation of the temperature field; such crossover (actu- 
ally, the oscillatory decay of the temperature field pertur- 
bations) is inconsistent with parabolic Fourier equation. 
So, if the size of the system is close to this critical scale, 
more exact macroscopic equations should be used for de- 
scription of the non-stationary heat conduction. In both 
models studied the critical size decreases with the tem- 
perature increase. As for the CV equation itself, in the 
FPU chain this equation clearly contradicts the simula- 
tions for the short-wave perturbations of the temperature 
field. In the chain of rotators it seems to be inconsistent 
with the simulations in the vicinity of the critical wave- 
length, however is more or less justified for longer and 
shorter modes. One can speculate that this difference 
between two models is related to their difference with 
respect to the stationary heat conduction - saturating 
versus size dependent behavior of the heat conduction 
coefficient [9, 10, 11, 12, 13]. 
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